Testing Plotting Functions

TO DO: - add plot with eTIV corrections - add CI’s

Prep

Let’s load a couple random-ish gamlss models and the dataframe they’re built on

set.seed(12345)
cn_df <- fread("/Users/megardn/Desktop/BGD_Repos/puberty_braincharts/cubic_scripts/bfp_iter/control_filt_og_data.csv", stringsAsFactors = TRUE, na.strings = "")


GMV.int <- readRDS("/Users/megardn/Desktop/BGD_Repos/puberty_braincharts/cubic_scripts/bfp_iter/gamlss_objs/GMV/best_GMV.m3s0_D_mod.rds")

sGMV.re <- gamlss(formula = sGMV ~ bfp(log_age, powers = c(-2, -2, 3), shift = 0, scale = 1) + sex + fs_version + re(random = ~ 1|study) + sex:bfp(log_age, powers = c(1, 2, 2), shift = 0, scale = 1),
                              sigma.formula = ~ bfp(log_age, powers = c(-2, -2), shift = 0, scale = 1) + sex + re(random = ~ 1|study),
                              nu.formula = ~ 1, family = GG, data = na.omit(cn_df),
                              control = gamlss.control(n.cyc = 200), trace = FALSE)

## GAMLSS-RS iteration 1: Global Deviance = 1218118 
## GAMLSS-RS iteration 2: Global Deviance = 1218033 
## GAMLSS-RS iteration 3: Global Deviance = 1218022 
## GAMLSS-RS iteration 4: Global Deviance = 1218017 
## GAMLSS-RS iteration 5: Global Deviance = 1218017 
## GAMLSS-RS iteration 6: Global Deviance = 1218017 
## GAMLSS-RS iteration 7: Global Deviance = 1218017 
## GAMLSS-RS iteration 8: Global Deviance = 1218017 
## GAMLSS-RS iteration 9: Global Deviance = 1218017 
## GAMLSS-RS iteration 10: Global Deviance = 1218017 
## GAMLSS-RS iteration 11: Global Deviance = 1218017 
## GAMLSS-RS iteration 12: Global Deviance = 1218017 
## GAMLSS-RS iteration 13: Global Deviance = 1218017

Call the script holding the plot functions we want to test:

source("plotting_functions.R")

Test

Set-up Functions

These functions are called from the plotting functions.

sim.data() - takes the dataframe you built your GAMLSS model on and creates a 2 new dfs with simulated data (male vs female participants) across the age-range, assigning fs_version and study values to whatever is most common in the original df. Expects input df to have log_age, fs_version, and study. Also preps x-axis labels and defines which centiles you’ll be plotting.

Returns object sim, which is a list of objects.

sim <- sim.data(cn_df)
names(sim)

## [1] "ageRange"            "dataToPredictM"      "dataToPredictF"     
## [4] "tickMarks_log"       "tickLabels_log"      "tickMarks_unscaled" 
## [7] "tickLabels_unscaled" "desiredCentiles"

centile_predict() - predicts centiles based on df simulated by sim.data() using the predictAll() and qGG() functions. Calculates centiles, 50th centile peak values, and age at peaks separately on male and female dfs and returns each, as well as an averaged effect across sexes.

Takes GAMLSS model obj and objects returned from sim.data(). Returns object pred, which is a list of objects.

pred <- centile_predict(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)

## new prediction 
## new prediction 
## new prediction 
## new prediction

names(pred)

##  [1] "fanCentiles"   "fanCentiles_M" "fanCentiles_F" "peak"         
##  [5] "peak_age"      "peak_M"        "peak_age_M"    "peak_F"       
##  [9] "peak_age_F"    "M_mu"          "M_sigma"       "M_nu"         
## [13] "F_mu"          "F_sigma"       "F_nu"

centile_predict.corrected() - same as centile_predict but corrects out the estimated effects of fs_version (from mu term) and study (from mu & sigma terms) parameters.

Expects fs_version to be a fixed effect and study to be a random effect fit using re() function!!!

pred.cor <- centile_predict.corrected(sGMV.re, sim$dataToPredictM, sim$dataToPredictF, sim$ageRange, sim$desiredCentiles)

## new prediction 
## new prediction 
## new prediction 
## new prediction

names(pred.cor)

##  [1] "fanCentiles"   "fanCentiles_M" "fanCentiles_F" "peak"         
##  [5] "peak_age"      "peak_M"        "peak_age_M"    "peak_F"       
##  [9] "peak_age_F"    "M_mu"          "M_sigma"       "M_nu"         
## [13] "F_mu"          "F_sigma"       "F_nu"

GGalt.variance - copied from Simon’s Nature paper repo

var <- GGalt.variance(pred$M_mu, pred$M_sigma, pred$M_nu)

Plotting

makeCentileFan() - basic centile fan plotting that averages across sex and predicts on Mode(fs_version) and Mode(study) of original data the gamlss was modeled on. Expects GAMLSS model, phenotype being modeled, and the name of the original df.

age_transformed parameter set to TRUE or FALSE

makeCentileFan(GMV.int, "GMV", cn_df, TRUE, "sex")

## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)

makeCentileFan(GMV.int, "GMV", cn_df, FALSE, "sex")

## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)

makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan(sGMV.re, "sGMV", cn_df, FALSE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan_sex_overlay() - same as makeCentileFan but with separate centile lines for males and females

makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, FALSE, "sex")

## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)

makeCentileFan_sex_overlay(GMV.int, "GMV", cn_df, TRUE, "sex")

## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)

makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, FALSE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan_sex_overlay(sGMV.re, "sGMV", cn_df, TRUE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan.corrected() - centile fan plot that averages across sexes, correcting for fs_version and study effects in both centiles and data points plotted.

Expects fs_version to be a fixed effect and study to be a random effect fit using re() function!!!

age_transformed parameter set to TRUE or FALSE

makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

#makeCentileFan(sGMV.re, "sGMV", cn_df, TRUE, "sex")
makeCentileFan.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan_sex_overlay.corrected() - same as makeCentileFan.corrected but with separate centile lines for males and females

makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, TRUE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

makeCentileFan_sex_overlay.corrected(sGMV.re, "sGMV", cn_df, FALSE, "sex")

## new prediction 
## new prediction 
## new prediction 
## new prediction

plot.gamlss.var() - plots variance from predicted GAMLSS model separately for males and females.

age_transformed parameter set to TRUE or FALSE

plot.gamlss.var(sGMV.re, "sGMV", cn_df, TRUE)

## new prediction 
## new prediction 
## new prediction 
## new prediction

plot.gamlss.var(GMV.int, "GMV", cn_df, FALSE)

## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)